#include "field.h"

void field_add( Field *f_a, const Field *f_b )
{
    CHECK_FIELD(f_a);
    CHECK_FIELD(f_b);
    CHECK_FIELDS(f_a,f_b);
    {
        double       **a = ARRAY_OF(f_a);
        const double **b = CONST_ARRAY_OF(f_b);
        const int      nx = f_a->nx;
        const int      ny = f_a->ny;
        int i,j;
        
        for(j=0;j<ny;++j)
        {
            for(i=0;i<nx;++i)
            {
                a[j][i] += b[j][i];
            }
        }
    }
}

void field_sub( Field *f_a, const Field *f_b, const double lam )
{
    CHECK_FIELD(f_a);
    CHECK_FIELD(f_b);
    CHECK_FIELDS(f_a,f_b);
    {
        double       **a = ARRAY_OF(f_a);
        const double **b = CONST_ARRAY_OF(f_b);
        const int      nx = f_a->nx;
        const int      ny = f_a->ny;
        int i,j;
        
        for(j=0;j<ny;++j)
        {
            for(i=0;i<nx;++i)
            {
                //const double tmp = a[j][i];
                a[j][i] -= lam*b[j][i];
                /*
                if( a[j][i] <= 0 )
                {
                    printf("rho=%g / drho=%g / dt=%g\n", a[j][i], b[j][i], lam);
                    exit(1);
                }*/
            }
        }
    }
}

void field_muladd( Field *f_a, const double lam, const double eta, const Field *f_b, const Field *f_c, const Field *f_d, const Field *f_dvg)
{
    CHECK_FIELD(f_a);
    CHECK_FIELD(f_b);
	CHECK_FIELD(f_c); 
    CHECK_FIELDS(f_a,f_b);
	CHECK_FIELDS(f_a, f_c);
    {
        double       **a   = ARRAY_OF(f_a);
        const double **b   = CONST_ARRAY_OF(f_b);
		const double **c   = CONST_ARRAY_OF(f_c);
		const double **d   = CONST_ARRAY_OF(f_d);
		const double **dvg = CONST_ARRAY_OF(f_dvg); 
        const int      nx  = f_a->nx;
        const int      ny  = f_a->ny;
        int i,j;
        
        for(j=0;j<ny;++j)
        {
            for(i=0;i<nx;++i)
            {
                a[j][i] += lam * (eta*b[j][i]-c[j][i]-d[j][i]+dvg[j][i]);
            }
        }
    }
}

void field_muladd_Ju(Field *f_a, const double lam, const Field *f_b)
{
	CHECK_FIELD(f_a);
	CHECK_FIELD(f_b);
	CHECK_FIELDS(f_a, f_b);
	{
		double **a = ARRAY_OF(f_a);
		const double **b = CONST_ARRAY_OF(f_b);
		const int nx = f_a->nx; 
		const int ny = f_a->ny; 
		int i,j; 
		
		for(j=0;j<ny;++j)
        {
            for(i=0;i<nx;++i)
            {
                a[j][i] += lam * b[j][i];
            }
        }
	}
}

